System and method for isogeometric analysis

ABSTRACT

A system and method for carrying out an isogeometric analysis process includes accessing the CAD model of the object from the memory. The process also includes analyzing the CAD model to generate a pre-refinement geometric map of the CAD object that has a smoothness projected to maintain a consistency of a geometric map based on the pre-refinement geometric map during a refinement of the mesh. The process further includes refining the mesh based on the pre-refinement geometric map to converge toward a refinement criteria associated with the CAD model.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

This invention was made with government support under 1435072 awarded bythe National Science Foundation. The government has certain rights inthis invention.

REFERENCE TO RELATED APPLICATION

N/A.

BACKGROUND OF THE INVENTION

The present disclosure relates, generally, to computer aided designsystems and methods.

Finite element analysis (FEA) is a computational tool commonly used byengineers in designing parts and systems using computer aided design(CAD) software. FEA allows testing of mechanical properties of acomputer designed part so that the part can be modified if it doesn'tmeet the required mechanical specifications. One of the limitations ofFEA is that it requires use of an approximated version of thecomputer-modeled part. The approximated geometry tends to be blockierthan the actual design, but this is required in order for conventionalFEA to break the part into tiny elements, which are each analyzed.Curved surfaces do not easily break up into the types of elements neededfor modern FEA. As a result the results are approximations of theperformance of the actual part, which in some cases are pretty close,but which in other cases can be quite far off. These elements areconnected together in the form of a mesh.

Isogeometric analysis is an analysis approach introduced by Hughes etal. in T. Hughes, J. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD,finite elements, NURBS, exact geometry and mesh refinement, Computermethods in applied mechanics and engineering 194 (39) (2005) 4135-4195,which is incorporated herein by reference in its entirety. Inisogeometric analysis, the same basis functions used to representgeometric models, such as Non-Uniform Rational B-Splines (NURBS), arealso used to approximate field variables in solving partial differentialequations (PDEs). Due to the same basis used in geometric representationand in solution approximation, isogeometric analysis eliminates thegeometric approximation error commonly occurred in classical FEAprocedures. Once an initial mesh is constructed, refinements can also beimplemented and an exact geometry is maintained at all levels withoutthe necessity of interaction with the CAD system. Another advantage ofisogeometric analysis is its computational efficiency on a per-nodebasis over classical C⁰ Lagrange polynomial based finite element. Thehigher continuity of the NURBS basis has been demonstrated tosignificantly improve the numerical efficiency and accuracy on a pernode basis in many areas including structural analysis, fluidsimulation, and shape optimization.

Isogeometric analysis techniques relying on a basis other than NURBShave also been developed. To overcome the limitation of the tensorproduct structure of NURBS in local mesh refinement, methods based onsubdivision solids and T-splines have been developed recently and havebeen successfully used in isogeometric analysis. The introduction ofT-junction in T-splines allows T-splines to represent complex shapes ina single patch and permit local refinement. On the other hand,challenges exist with respect to analysis-suitable T-splines, such ashow to obtain efficient local refinement and effective treatment ofso-called extraordinary points.

Recently, triangular Bézier splines have emerged as a powerfulalternative to shape modeling and isogeometric analysis due to theirflexibility in representing shapes of complex topology and their higherorder of continuity. Local refinement can also be implemented withoutany great difficulty. Various sets of basis functions have been definedon triangulations using bivariate spline functions. Some bivariatesplines have also been effectively applied in solving PDEs, includingsome where quadratic Powell-Sabin (PS) splines with C¹ smoothness areconsidered. A locally-supported basis is constructed by normalizing thepiecewise quadratic PS B-splines and it is cast in the Bernstein-Bézierform. Such basis has global C¹ smoothness and is used to approximate thesolution of PDEs. More generalized C^(r) basis and elements includingPS, Clough-Tocher (CT), and polynomial macro-elements based on rationaltriangular Bézier splines (rTBS) have also been introduced and appliedsuccessfully in isogeometric analysis on triangulations.

Referring to FIG. 1, cubic C¹ smooth basis functions with CTmacro-elements are illustrated within the context of a rTBS basedisogeometric analysis. The given physical domain is triangulated into aset of C¹ smooth Bézier elements, which are mapped from the parametricmesh. More particularly, free control points and domain points 10 aredetermined by dependent control points and domain points 12,respectively under the continuity constraints. Thus, the C¹ basisfunctions ψ are constructed as linear combinations of the C⁰ Bernsteinbasis ϕ, under the continuity constraints. The resulting analysis hasshown to be efficient, accurate, and convergent. However, optimalconvergence in h-refinement has only been achieved for C⁰ elements andthe convergence rate is sub-optimal for C^(r) elements.

Thus, it would be desirable to have a system and method for creatingmeshes from models with a controlled amount of and, preferably, withoutany approximation. Furthermore, it would be desirable that such systemand method demonstrate effective and efficient convergence in theanalysis.

SUMMARY OF THE INVENTION

In accordance with one aspect of the present disclosure, a system isprovided for creating a mesh from a computer aided design model duringan isogeometric analysis process. The system includes a memory havingaccess to a computer aided design (CAD) model of an object. The systemalso includes a processor configured to carry out an isogeometricanalysis process. The processor is configured to access the CAD model ofthe object from the memory. The processor is also configured to analyzethe CAD model to generate a pre-refinement geometric map of the CADobject that has a smoothness projected to maintain a consistency of amesh based on the pre-refinement geometric map during a refinement ofthe mesh. The processor is further configured to refine the mesh basedon the pre-refinement geometric map to converge toward a refinementcriteria associated with the CAD model.

In accordance with another aspect of the present invention, a method isprovided for creating a mesh from a computer aided design model duringan isogeometric analysis process. The method includes accessing to acomputer aided design (CAD) model of an object. The method also includesanalyzing the CAD model to generate a pre-refinement geometric map toserve as a map between a parametric mesh and a physical mesh of thecomputer aided design model. To do so, the method includes generatingthe pre-refinement geometric map to have smoothness projected tomaintain consistency of the geometric map between the parametric mesh orthe physical mesh during refinement by determining control points thatdo not need to be relocated during refinement of the parametric mesh orthe physical mesh because the control points satisfy a continuityconstraint. The method further includes refining the parametric mesh orthe physical mesh based on the pre-refinement geometric map to convergetoward a refinement criteria associated with the CAD model.

Additional features and advantages of the present invention will beapparent from the following detailed description taken in conjunctionwith the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic illustration of a rational triangular Béziersplines (rTBS) based isogeometric analysis.

FIG. 2A is a schematic illustration of domain points of the Bézierordinates b_(ijk) in {v₁, v₂, v₃}.

FIG. 2B is a graphic illustration of triangular Bézier patch b(ξ).

FIG. 3A is a graphic illustration of triangular Bézier patches with C¹continuity constraints, specifically, two domain triangles with C¹constraints on Bézier ordinates.

FIG. 3B is a graphic illustration of triangular Bézier patches with C¹continuity constraints, specifically, two Bézier patches with C¹continuity.

FIG. 4A is a graphic illustration of a cubic C¹ mesh using Clough-Tocher(CT) split.

FIG. 4B is a graphic illustration of a quadratic C¹ mesh usingPowell-Sabin (PS) split.

FIG. 5 is a flow chart setting forth steps of a process for automaticdomain parameterization.

FIG. 6A is a graphic illustration of an input domain with NURBSboundary.

FIG. 6B is a graphic illustration of the input with Bézier extraction.

FIG. 6C is a graphic illustration of the input with polygonal parametricdomain.

FIG. 6D is a graphic illustration of the input with parametric mesh withdomain points.

FIG. 6E is a graphic illustration of the input with physical mesh afterboundary replacement.

FIG. 6F is a graphic illustration of the input with local subdivision onthe circle.

FIG. 6G is a graphic illustration of the input with a new parametricdomain.

FIG. 6H is a graphic illustration of the input with a new parametricmesh.

FIG. 6I is a graphic illustration of the input with a new physical meshafter boundary replacement.

FIG. 6J is a detailed view of boundary replacement for a portion of FIG.6(e)

FIG. 6K is a detailed view of boundary replacement for a portion of FIG.6(i)

FIG. 7 is a graphic illustration showing Πuƒ is the push-forward of therTBS projector Π_(S)(ƒ∘G), where ƒϵL²(Ω) and ƒ∘GϵL²({circumflex over(β)}).

FIG. 8A is a graphic illustration of an initial quadratic parametricmesh {circumflex over (T)}⁰.

FIG. 8 B is a graphic illustration of an initial quadratic physical meshT⁰.

FIG. 9A is a graphic illustration of a refinement of the elements ofFIG. 8B without imposing the continuity constraints illustrating T_(ct*)⁰.

FIG. 9B is a graphic illustration of a refinement of the elements ofFIG. 8B without imposing the continuity constraints illustratingT_(u,ct*) ⁰.

FIG. 9C is a graphic illustration of a refinement of the elements ofFIG. 8B without imposing the continuity constraints illustratingT_(u,u,ct*) ⁰.

FIG. 9D is a graphic illustration of a refinement of the elements ofFIG. 8B with C¹ continuity constraints illustrating T_(ct) ¹.

FIG. 9E is a graphic illustration of a refinement of the elements ofFIG. 8B with C¹ continuity constraints illustrating T_(u,ct) ¹.

FIG. 9F is a graphic illustration of a refinement of the elements ofFIG. 8B with C¹ continuity constraints illustrating T_(u,u,ct) ¹.

FIG. 10A is a graphic illustration of an initial mesh, T^(1,2).

FIG. 10B is a graphic illustration of the initial mesh of FIG. 10A withone refinement, T_(u) ^(1,2).

FIG. 10C is a graphic illustration of the refined mesh of FIG. 10A withtwo refinements, T_(u,u) ^(1,2) (what do the colors of the points mean?)

FIG. 11A is a graphic illustration showing a cubic C¹ mesh obtained byDC with CT macro-elements, where the space s₃ ¹ is well defined on acubic mesh with CT macro-elements.

FIG. 11B is a graphic illustration showing a C¹ mesh after uniformlyrefining mesh of FIG. 11A, where the dimension of the space s₃ ¹ on thistriangulation is unknown.

FIG. 11C is a graphic illustration showing a CT split is performed onthe mesh of FIG. 11B to construct a set of stable local basis thatdefines s₃ ¹.

FIG. 12A is a graphic illustration showing a problem domain.

FIG. 12B is a graphic illustration showing an initial parametric meshrelative to the problem domain of FIG. 12A.

FIG. 12C is a graphic illustration showing an initial physical meshrelative to the problem domain of FIG. 12A.

FIG. 12D is a graphic illustration showing is a quadratic mesh, S₂ ⁰,relative to the problem domain of FIG. 12A.

FIG. 12E is a graphic illustration showing a cubic mesh, S₃ ⁰, relativeto the problem domain of FIG. 12A.

FIG. 12F is a graphic illustration showing a quintic mesh, S₅ ⁰,relative to the problem domain of FIG. 12A.

FIG. 13A is a graphic illustration showing a physical mesh relative tothe problem of FIG. 12A, S₂ ¹(T_(ps)) by DC, dim=36.

FIG. 13B is a graphic illustration showing a physical mesh relative tothe problem of FIG. 12A, S₂ ¹(T_(ps)) by GE, dim=36.

FIG. 13C is a graphic illustration showing a physical mesh relative tothe problem of FIG. 12A, S₂ ¹(T_(ps)) by DC, dim=60.

FIG. 13D is a graphic illustration showing a physical mesh relative tothe problem of FIG. 12A, S₃ ¹(T_(ct)) by GE, dim=60.

FIG. 13E is a graphic illustration showing a physical mesh relative tothe problem of FIG. 12A, S₅ ^(1,2)(T) by DC, dim=96.

FIG. 13F is a graphic illustration showing a physical mesh relative tothe problem of FIG. 12A, S₅ ¹(T) by GE, dim=120.

FIG. 13G is a graphic illustration showing a physical mesh relative tothe problem of FIG. 12A, S₅ ^(2,3)(T_(ps)) by DC, dim=120.

FIG. 13H is a graphic illustration showing a physical mesh relative tothe problem of FIG. 12A, S₅ ²(T_(ps)) by GE, dim=252.

FIG. 14A is a graph showing convergence rates in C⁰ space.

FIG. 14B is a graph showing convergence rates in C¹ space.

FIG. 14C is a graph showing convergence rates in C² space.

FIG. 15A is a graphic illustration of an L-shaped domain with threeholes.

FIG. 15B is a graphic illustration of an initial quadratic parametricmesh of the domain of FIG. 15A.

FIG. 15C is a graphic illustration of an initial quadratic physical meshof the domain of FIG. 15A.

FIG. 16A is a graphic illustration showing a parametric mesh in S₂¹(T_(ps)) with basis obtained by DC with PS macro-elements.

FIG. 16B is a graphic illustration showing a C¹ physical mesh obtainedby DC with PS macro-elements.

FIG. 16C is a graphic illustration showing a parametric mesh in S₃¹(T_(ct)) with basis obtained by DC with CT macro-elements.

FIG. 16D is a graphic illustration showing a C¹ physical mesh obtainedby DC with CT macro-elements.

FIG. 16E is a graphic illustration showing a parametric mesh in S₅ ¹(T)with basis obtained by GE.

FIG. 16F is a graphic illustration showing a C¹ physical mesh obtainedby GE.

FIG. 16G is a graphic illustration showing a parametric mesh in S₅^(1,2)(T) with basis obtained by DC.

FIG. 16H is a graphic illustration showing a C² physical mesh in S₅^(1,2)(T).

FIG. 17A is a graph showing error measured in the L²-norm vs. meshparameter from refined elements in FIGS. 16A-16H to illustrateconvergence rates in C⁰ space.

FIG. 17B is a graph showing error measured in the L²-norm vs. meshparameter from refined elements in FIGS. 16A-16H to illustrateconvergence rates in C¹ space with nested refinement sequences.

FIG. 17C is a graph showing Error measured in the L²-norm vs. meshparameter in C¹ spaces for refinement sequences with inconsistentgeometric map, establishing that the optimal convergence rates are notachieved.

FIG. 18 is a graphic illustration of an elastic plate with a circularhole, where L is the length of the edge, R is the radius of the circle,and τ is the thickness of the plate, and where E and ν represent theYoung's modulus and Poisson ratio respectively.

FIG. 19A is a graphic illustration of an initial parametric mesh of thephysical domain of FIG. 18.

FIG. 19B is a graphic illustration of an initial physical mesh of thephysical domain of FIG. 18.

FIG. 19C is a graphic illustration of a smoothed parametric mesh of thephysical domain of FIG. 18.

FIG. 20A is a graphic illustration of a parametric mesh in S₃ ⁰(T) ofthe physical domain of FIG. 18.

FIG. 20B is a graphic illustration of a physical mesh in S₃ ⁰(T) of thephysical domain of FIG. 18.

FIG. 20C is a graphic illustration of a parametric mesh in S₅ ⁰(T) ofthe physical domain of FIG. 18.

FIG. 20D is a graphic illustration of a physical mesh in S₅ ⁰(T) of thephysical domain of FIG. 18.

FIG. 20E is a graphic illustration of a parametric mesh in S₂ ¹(T_(ps))with basis obtained by DC with PS macro-elements of the physical domainof FIG. 18.

FIG. 20F is a graphic illustration of a C¹ physical mesh obtained by DCwith PS macro-elements of the physical domain of FIG. 18.

FIG. 20G is a graphic illustration of a parametric mesh in S₃ ¹(T_(ct))with basis obtained by DC with CT macro-elements of the physical domainof FIG. 18.

FIG. 20H is a graphic illustration of a C¹ physical mesh obtained by DCwith CT macro-elements of the physical domain of FIG. 18.

FIG. 20I is a graphic illustration of a parametric mesh S₅ ^(1,2)(T)with basis obtained by DC of the physical domain of FIG. 18.

FIG. 20J is a graphic illustration of a C² physical mesh in S₅ ^(1,2)(T)obtained by GE of the physical domain of FIG. 18.

FIG. 20K is a graphic illustration of a parametric mesh in S₅ ¹(T) withbasis obtained by GE of the physical domain of FIG. 18.

FIG. 20L is a graphic illustration of a C¹ physical mesh in S₅ ¹(T)obtained by GE of the physical domain of FIG. 18.

FIG. 20M is a graphic illustration of a parametric mesh in S₅^(2,3)(T_(ps)) with basis obtained by DC with PS macro-elements of thephysical domain of FIG. 18.

FIG. 20N is a graphic illustration of a C³ physical mesh in S₅^(2,3)(T^(ps)) obtained by GE of the physical domain of FIG. 18.

FIG. 20O is a graphic illustration of a parametric mesh in S₅ ²(T_(ps))with basis obtained by GE with PS macro-elements of the physical domainof FIG. 18.

FIG. 20P is a graphic illustration of a C² physical mesh in S₅ ²(T_(ps))obtained by GE of the physical domain of FIG. 18.

FIG. 21A is a graph showing convergence rates in C⁰ space.

FIG. 21B is a graph showing convergence rates in C¹ space.

FIG. 21C is a graph showing convergence rates in C² space S₅ ²(T_(ps))(GE) and superspline space S₅ ^(2,3)(T_(ps)) (DC).

FIG. 22 is a graph showing error measured in the L²-norm of stress vs.the number of nodes from refinements of three quintic elements, S₅ ⁰, S₅¹, S₅ ^(2,3), in FIGS. 20C/20D, FIG. 20K/20L and FIG. 20M/20N

FIG. 23 is a graph showing error measured in the L²-norm of stress vs.mesh parameter for refinement sequences with inconsistent geometric map,and illustrating that the convergence rates are remarkably lower than inFIGS. 21A-21C.

FIG. 24A is a graphic illustrating an initial S₃ ¹(T_(ct)) mesh.

FIG. 24B is a graphic illustrating a local refinement 1.

FIG. 24C is a graphic illustrating a local refinement 2.

FIG. 24D is a graphic illustrating a local refinement 3.

FIG. 25 is a graph showing a comparison of error per degree of freedombetween uniform and local refinement in space S₃ ¹(T_(ct)) for the platehole problem.

FIG. 26 is a block diagram of a non-limiting of an example of one systemfor implementing the systems and methods described herein.

FIG. 27 is a flowchart setting forth steps of one example of a processin accordance with the present disclosure.

DETAILED DESCRIPTION OF THE INVENTION

As will be described, the present disclosure provides an approach thatcan provide convergence for all C^(r) rTBS elements based isogeometricanalysis in the context of h-refinement, even relative to a criteria foroptimal convergence. Approximation power and convergence rates are usedto evaluate the performance of numerical schemes for solving PDEs. TheNURBS space has been proved to have full approximation power and deliveran optimal rate of convergence as the classical finite element spaces.Similar approximation estimates that are optimal with respect to thepolynomial degree of the underlying spline space have also beendeveloped for T-splines that are defined on a particular two-patchstructure. However, no convergence results have been reported forT-splines defined on more generalized T-meshes that includeextraordinary points, which usually require additional continuityconstraints on the surrounding control points to achieve G¹ continuity.

As will follow, an approach to rTBS based isogeometric analysis isdescribed. An initial coarse parametric mesh is first refined intoelements that are sufficiently small for analysis and the global C^(r)smooth basis is then constructed by imposing continuity constraints onadjacent triangles in the refined parametric mesh. Based on the C^(r)basis, the C^(r) geometric map is obtained. Although such arefine-then-smooth approach does lead to C^(r) stable basis that issufficient for analysis, the resulting geometric map may not beconsistent before and after the C^(r) constraints are imposed. Theinconsistency in geometric maps leads to deteriorated convergence rate.

The present disclosure recognizes that the reason for such inconsistencyin geometric map is as follows. The C^(r) basis is obtained viacontinuity constraints on domain points in the parametric mesh. With theC^(r) basis, some domain points are free and other domain points aredependent on these free points. For the geometric map that maps theparametric mesh to the physical mesh, the control points correspondingto the free domain points are chosen as free control points. If theremaining dependent control points do not satisfy the same continuityconstraints, they would have to be relocated to satisfy the constraintsto ensure the map is C^(r). Such relocation of dependent control pointsleads to a change of the geometric map.

To overcome such inconsistency in the geometric map in h-refinement withthe refine-then-smooth approach, a three-step approach is provided. Aswill described, this approach can to achieve convergence in IGA withC^(r) rTBS elements. In one non-limiting example, this approach canachieve convergence with respect to objective, optimal criteria forconvergence. In particular, a pre-refinement geometric map isconstructed that possesses sufficient smoothness to maintain theconsistency of the geometric map for subsequent refinements. From thepre-refinement smooth geometric map, the mesh is uniformly refined. Themacro-element techniques is used to obtain stable C^(r) smooth basis foranalysis. In this “smooth-refine-smooth” approach, the smoothness in thefirst step is used for geometric reason, so that the resulting geometricmap stays consistent during the mesh refinement. The smoothness in thelast step provides a stable basis over triangulation for analysis. Theprovided smooth-refine-smooth approach can, as a non-limiting example,provide optimal convergence in h-refinement for all types of C^(r) rTBSelements. This approach is also applicable to some supersplines S₅^(r,ρ), (ρ>r), where some vertices or edges in macro-triangles possesshigher order C^(ρ) smoothness than the global C^(r) smoothness. In suchcases, the smoothness of the pre-refinement geometric map should beC^(ρ).

Bézier Triangles

NURBS has been widely used as a standard to represent curves andsurfaces in CAD systems. Each knot span of a NURBS curve corresponds toa Bézier curve that is defined through Bernstein basis functions. A d-thdegree Bernstein polynomial is defined as:

$\begin{matrix}{{{B_{{ij},d}(\xi)} = {\begin{pmatrix}d \\{i,j}\end{pmatrix}{\xi^{i}\left( {1 - \xi} \right)}^{d - i}}},{\xi \in \left\lbrack {0,1} \right\rbrack},} & {(1);}\end{matrix}$

where

${\left( \frac{d}{i,j} \right) = \frac{d!}{{i!}{j!}}},{{i + j} = {d.}}$Accordingly a d-th degree bivariate Bernstein polynomial is defined as:

$\begin{matrix}{{{B_{i,d}(\xi)} = {\frac{d!}{{i!}{j!}{k!}}\gamma_{1}^{i}\gamma_{2}^{j}\gamma_{3}^{k}}},\mspace{14mu}{{i} = {{i + j + k} = d}},} & {(2);}\end{matrix}$

where i represents a triple index (i,j,k) and (γ₁,γ₂,γ₃) is thebarycentric coordinate of a point ξϵ

. Every point ξ=(ξ₁,ξ₂) in a fixed triangle with vertices v₁, v₂, v₃ϵ

² can be written uniquely in the form:ξ=γ₁ v ₁+γ₂ v ₂+γ₃ v ₃,  (3);

with γ₁+γ₂+γ₃=1.

It has been shown that the set

{B_(i, d)}_(i = d)is a basis for the space of degree d bivariate polynomials

_(d). A triangular Bézier patch can be defined as:

$\begin{matrix}{{{b(\xi)} = {\sum\limits_{{i} = d}{p_{i}{B_{i,d}(\xi)}}}},} & {(4);}\end{matrix}$

where p_(i) represents a triangular array of control points. A rationalBézier triangle can be defined similarly as:

$\begin{matrix}{{{b(\xi)} = {\sum\limits_{{i} = d}{p_{i}{\phi_{i,d}(\xi)}}}},} & {(5);}\end{matrix}$

with ϕ_(i,d) being the rational Bernstein basis:

$\begin{matrix}{{\phi_{i,d} = {\frac{w_{i}B_{i,d}}{\sum\limits_{{i} = d}{w_{i}B_{i,d}}} = \frac{w_{i}B_{i,d}}{w}}},} & {(6);}\end{matrix}$

where w_(i) are the weights associated with the control points p_(i).

Bivariate Bernstein polynomials can be used to define a polynomialfunction f of degree d over a triangle τ={v₁, v₂, v₃} as:

$\begin{matrix}{{{f(\xi)} = {\sum\limits_{{i} = d}{b_{i}{\phi_{i,d}(\xi)}}}},} & (7)\end{matrix}$

where ϕ_(i,d) are the rational Bernstein basis polynomials associatedwith τ. The b_(i) are called the Bézier ordinates of f and theirassociated set of domain points is defined as:

$\begin{matrix}{D_{d,T} = {\left\{ {{q_{ijk} = \frac{{iv}_{1} + {jv}_{2} + {kv}_{3}}{d}},{{i + j + l} = d}} \right\}.}} & (8)\end{matrix}$

For example, FIGS. 2A and 2B illustrate an example of the associateddomain points of the Bézier ordinates (2A) and triangular Bézier patch(2B). Two polynomials f and {hacek over (ƒ)} of degree d join r timesdifferentially across the common edge of two triangles τ={v₁, v₂, v₃}and {tilde over (τ)}={v₄, v₃, v₂}, if and only if:

$\begin{matrix}{{{{\overset{\sim}{b}}_{\rho,j,k} - {\sum\limits_{{u + v + \kappa} = \rho}{\frac{\rho!}{{\mu!}{v!}{\kappa!}}b_{\mu,{k + v},{j + \kappa}}\gamma_{1}^{\mu}\gamma_{2}^{v}\gamma_{3}^{\kappa}}}} = 0},{{j + k + \rho} = d},{\rho = 0},\ldots\mspace{14mu},r,} & {(9);}\end{matrix}$

where γ₁,γ₂,γ₃ are the barycentric coordinates of vertex v₄ with respectto triangle τ. For example, FIGS. 3A and 3B illustrate two triangularBézier patches with C¹ continuity constraints. In FIG. 3A, two domaintriangles are illustrated with C¹ constraints on Bézier ordinates and,in FIG. 3B, two Bézier patches are illustrated with C¹ continuity. Inboth, free nodes 30, whose values can be freely chosen, are shown. Also,dependent nodes 32 are determined by the free nodes 30 through thecontinuity constraints. The triangles where continuity constraints areimposed 34 are further illustrated. As can be seen in FIG. 3B, thecontrol points in each shaded triangle pair 34 are coplanar. For bettervisualization of the underlying C^(d) patch, the control net in FIG. 3Bis shifted up slightly.

Splines on Triangulations

Consider a parametric domain {circumflex over (Ω)} and its triangulation{circumflex over (T)}. Then, introduce the spline spaces of piecewisepolynomials of degree d over {circumflex over (T)}:S _(d) ^(r)({circumflex over (T)})={(ƒϵC ^(r)({circumflex over(Ω)}):ƒ|_(τ) ϵP _(d)∀_(τ) ϵ{circumflex over (T)}},  (10)

where τ is an arbitrary triangle in {circumflex over (T)} and r is thecontinuity order of the spline over {circumflex over (Ω)}. In addition,if the spline has higher smoothness at some vertices or across someedges, the spline can be considered a superspline and the associatedspace denoted as:S _(d) ^(r,ρ)({circumflex over (T)})={ƒϵS _(d) ^(r)({circumflex over(T)}):ƒϵC ^(ρ) ^(v) (v)∀vϵV&ƒϵC ^(ρ) ^(e) (e)∀eϵE}  (11)

where V and E are the set of all vertices and edges respectively in{circumflex over (T)} and ρ:={ρ_(v)}_(vϵV)∪{ρ_(e)}_(eϵE) withr≤ρ_(v),ρ_(e)≤d for each vϵV and eϵE.

There are several approaches to obtain C^(r) spline spaces on atriangulated domain {circumflex over (Ω)}({circumflex over (T)}). Withrespect to the spaces S_(d) ^(r) and S_(d) ^(r,ρ) with fullapproximation power of d-th degree polynomials, condition equation (9)can be applied directly on the triangles, which requires the degree ofthe polynomial much higher than r, such as d≥3r+2. An alternativestrategy is to split each triangle in T into several micro-trianglesbefore imposing the continuity constraints on the micro-triangles. Theoriginal triangles are then called macro-triangles. These include theClough-Tocher (CT) split with polynomials of degree d≥3r for continuityr-odd and d≥3r+1 for r-even, and the Powell-Sabin (PS) split withpolynomials of degree

$d \geq \frac{{9r} - 1}{4}$for r-odd and of degree

$d \geq \frac{{9r} + 4}{4}$for r-even. For example, CT split may be used to obtain S₃ ¹ splinespace with cubic polynomials, and PS split may be used to obtain S₂ ¹,S₅ ² and S₅ ^(2,3) spline spaces with quadratic and quintic polynomials,respectively. The so-called polynomial macro-element technique may beused to obtain S₅ ¹ and S₅ ^(1,2) spline spaces with quintic polynomialswithout using any split technique.

Specifically, FIGS. 4A and 4B show the CT and PS splits, respectively,with corresponding free points 40 and dependent domain points 42. In theCT split illustrated in FIG. 4A, each vertex of a triangle in T isconnected with its centroid point to form three micro-triangles. Thisresulting triangulation is denoted as {circumflex over (T)}_(ct). In thePS split, each triangle is connected along lines from its incenter toeach of the three vertices and connected from incenter to a common edgebetween triangles. In addition, the middle of each boundary edge to theincenter of the associated triangle is connected, to thereby create sixmicro-triangles. In this implementation, the centroid point, instead ofthe incenter of each triangle, is used as the interior split point andthe resulting triangulation is denoted as {circumflex over (T)}_(ps).

Uniform refinement can also be performed as needed. For example, eachtriangle can be subdivided into four sub-triangles by connecting themiddle points of the edges. This kind of 1-to-4 split based uniformrefinement is used in our subsequent analysis of convergence during meshrefinement.

Automatic Domain Parametrization with C^(r) rTBS

As will be described, a method is provided for discretizing a physicaldomain Ω into a collection of rational Bézier triangles without anygeometric approximation error. Given an arbitrary 2D domain Ω and itsNURBS-represented boundary Γ, we seek a geometric map G(ξ),ξϵ{circumflex over (Ω)} such that the physical domain Ω is the image ofthe geometric map G(ξ) over a parametric domain {circumflex over (Ω)}where the physical boundary is reproduced by the map. In addition, thegeometric map G(ξ) is continuous and differentiable up to any desireddegree of continuity C^(r).

This can be achieved in three general steps that form a process 50, asillustrated with respect to FIG. 5. First, form a polygonal parametricdomain {circumflex over (Ω)} and its triangulation {circumflex over (T)}from the given physical domain Ω and the triangulation T₀, asillustrated at process block 51. Establish a C⁰ geometric map G₀ betweenthem: G₀:{circumflex over (Ω)}_(T)

Ω_(T) ₀ . Second, at process block 52, construct a set of C^(r) basisψ(ξ) on {circumflex over (Ω)}. Third, at process block 53, construct aC^(r) continuous triangulation T of Ω and establish a globally C^(r)geometric map G(ξ) from the parametric domain {circumflex over(Ω)}_({circumflex over (T)}) to the physical domain Ω_(T). As will bedescribed, this process can be improved with additional control againstpossible self-intersection in the physical mesh.

Construction of a C⁰ Geometric Map G₀

Referring to FIGS. 5 and 6A-6I, the process for constructing aparametric and physical mesh {circumflex over(Ω)}_({circumflex over (T)}) and Ω_(T) for a given input domain boundedby NURBS curves can be described. In FIG. 6, end points 60 and interiorcontrol points 62 of NURBS curves are shown along the boundaries, as arethe domain points 64 and control points of the mesh. Referring to FIGS.5 and 6A, at process block 54, given a domain Ω with NURBS boundarycurves of degree d, each NURBS curve can be subdivided into a set ofBézier curves via knot insertions, such as illustrated in FIG. 6B. Atprocess block 55, the end points of these Bézier curves are connected toform a polygonal parametric domain {circumflex over (Ω)}, as illustratedin FIG. 6C, and the domain {circumflex over (Ω)} is triangulated usingDelaunay triangulation to obtain {circumflex over(Ω)}_({circumflex over (T)}) and the associated domain points aregenerated according to Eq. 8, as illustrated in FIG. 6D. The quality ofthe parametric mesh may also be improved by using appropriate techniquessuch as a smoothing method, such as described in N. Jaxon, X. Qian,Isogeometric analysis on triangulations, Computer-Aided Design 46 (2014)45-57, which is incorporated herein by reference in its entirety. Atprocess block 56, the boundary control points of {circumflex over(Ω)}_({circumflex over (T)}) with corresponding control points of theBézier curves in the physical domain can be replaced to obtain atriangulation T₀ on the physical domain, such as illustrated in FIG. 6E.Thereafter, the physical mesh T₀ may be checked to see if there isself-folding. As illustrated in FIG. 6J, for example, this may happenbecause of the over recessed control points, the control polygon of thecurved boundary intersects with the other two boundaries of the element,causing self-folding of the mesh. If so, the Bézier boundary curve canbe subdivided where self-folding occurs and repeat steps (2)-(3) until avalid physical mesh T₀ is obtained, as illustrated in FIGS. 6F, 6G, and6H. On the other hand, as illustrated in 6K, the control polygon of thecurved boundary may be contained between the other two boundaries of theelement, the mesh is clear of self-folding and a valid physical mesh T₀is obtained, as illustrated in FIG. 6I and process block 58.

Construction of C^(r) spline basis ψ(ξ) on {circumflex over (T)}

Let

_(d,{circumflex over (T)})=denote the set of domain points fortriangulation {circumflex over (T)} as defined in equation (8), vϵ

_(d,{circumflex over (T)}) a domain point and b_(v) its ordinate. Apiecewise polynomial function ƒ(ξ)ϵS_(d), ξϵ

² can be expressed in terms of the rational C⁰ Bernstein basis ϕ andcorresponding nodal ordinates b

_(d,{circumflex over (T)}) as:

$\begin{matrix}{{f(\xi)} = {{\sum\limits_{i}{b_{i}{\phi_{i}(\xi)}}} = {b_{\mathcal{D}_{d,\hat{T}}}^{T}{{\phi(\xi)}.}}}} & (12)\end{matrix}$

Further, for a C^(r) continuous spline ƒϵS_(d) ^(r), the C^(r)smoothness conditions in equation (9) among the Bézier ordinates implythat we cannot assign arbitrary values to every coefficient of f.Instead, only certain coefficients corresponding to a reduceddetermining set of domain points

_(d,{circumflex over (T)})⊂

_(d,{circumflex over (T)}) can be assigned, and all remainingcoefficients may be determined by the smoothness conditions. When

_(d,{circumflex over (T)}) is the smallest set among all possibledetermining sets, we call

_(d,{circumflex over (T)}) a minimal determining set (MDS) and thedomain points in it free nodes. We define a set of basis ψ(ξ) for thespline space S_(d) ^(r)({circumflex over (T)}) in terms of these freenodes as:ψ(ξ)={ψ_(v) ϵS _(d) ^(r)({circumflex over (T)}):b _(v) b _(u)=δ_(v,u),∀u,vϵ

_(d,{circumflex over (T)})},  (13);

where ψ_(v) is the basis function at domain point v and δ_(v,u) is theKronecker delta.

The construction of such explicit MDS and, hence, the basis of theunderlying spline space is not a trivial task. One approach is directionconstruction through macro-elements. That is, according to theconnectivity of the triangle elements in a specific triangulation, onecan directly choose a set of free domain points based on which all otherdomain points are determined through necessary continuity constraints.This will be referred hereinto as the direct construction (DC) methodand the resulting spaces as macro-element spaces. Some examples includequintic C¹ polynomial macro-element space S₅ ^(1,2)({circumflex over(T)}), quadratic C¹ PS macro-element space S₂ ¹({circumflex over(T)}_(ps)), cubic C¹ CT macro-element space S₃ ¹({circumflex over(T)}_(ct)), and quintic C² PS macro-element space S₅ ^(2,3)({circumflexover (T)}_(ps)), where S₅ ^(1,2)({circumflex over (T)}) and S₅^(2,3)({circumflex over (T)}_(ps)) are in fact superspline spaces withS₅ ^(1,2)({circumflex over (T)}) having C² super-smoothness at everyvertex and S₅ ^(2,3)({circumflex over (T)}_(ps)) having C³super-smoothness at every vertex and splitting point of themacro-elements and across the three interior edges not connecting to thevertices of each macro-element.

Alternatively, the MDS can also be constructed by analyzing ahomogeneous linear system of the smoothness conditions using equation(9) for all pairs of triangles sharing an interior edge, which is:Ab

_(d,{circumflex over (T)}) =0  (14);

where A is a coefficient matrix depending on the geometry of the domaintriangles and b

_(d,{circumflex over (T)}) are n Bézier ordinates for the domain pointsin

_(d,{circumflex over (T)}). The dimension of the space S_(d)^(r)({circumflex over (T)}) is:dim S _(d) ^(r)({circumflex over (T)})=dim S _(d) ⁰({circumflex over(T)})−rank A.  (15).

Revealing the rank of A requires Gaussian elimination (GE). However,this is challenging in floating point arithmetic for geometry withdegeneracies, which would lead to increased rank deficiencies of A andcan be easily obscured by inexact computations. To overcome such issue,a modified Gaussian elimination procedure based on residual arithmetichas been presented. This method makes use of the fact that the verticesof the triangulation are pixels and the coordinates of the pixels areintegers.

In accordance with the present disclosure, besides the DC method, thestandard GE method can be used with complete pivoting to construct theMDS. For example, the GE method can be used in two types of situations.The first situation is for triangulations constructed throughmacro-elements, where the GE method can be used to identify stable basisfor analysis. Although it is not based on residual arithmetic, it workswell for the examples, such as will be provided. The second situation isfor general triangulations where the goal is to obtain smoothpre-refinement geometric map, rather than stable basis.

With the known MDS, after some manipulations on matrix A, equation (14)can be transformed to the form:

$\begin{matrix}{{b_{\mathcal{D}_{d,\hat{T}}} = {C^{T}b_{\mathcal{M}_{d,\hat{T}}}}};} & (16)\end{matrix}$

where C is called the continuity matrix. For the convenience of applyingDirichlet boundary conditions conveniently, a boundary MDS is alsoenforced, which means the complete boundary will be uniquely determinedby the free nodes on the boundaries only. This is accomplished byexchanging any free nodes with influence on the boundary with aconstrained boundary node. A free node appears in C as a column with asingle 1 that is otherwise all zeros. If a constrained boundary node isdependent on a free internal node, then by scaling this free basis row,and adding multiples of it to zero the boundary node's column, wereplace the internal free node by one on the boundary. Combiningequations (16) and (12), the C^(r) continuous function f now can beexpressed in terms of the free nodal ordinates

b_(ℳ_(d, T̂))as

$\begin{matrix}{{{f(\xi)} = {{b_{\mathcal{D}_{d,\hat{T}}}^{T}{\phi(\xi)}} = {{b_{\mathcal{M}_{d,\hat{T}}}^{T}C\;{\phi(\xi)}} = {b_{\mathcal{M}_{d,\hat{T}}}^{T}{\psi(\xi)}}}}};} & (17) \\\text{where:} & \; \\{{{\psi(\xi)} = {C\;{\phi(\xi)}}};} & (18)\end{matrix}$

is a set of global C^(r) basis functions composed as the linearcombinations of the C⁰ Bernstein basis ϕ(ξ).

Construction of C^(r) Geometric Map G(ξ)

After identifying the free control points p_(i) ^(ƒ) corresponding tothe free domain points in

_(d,{circumflex over (T)}), all the control points p for rTBS elementsin the physical mesh T may be overridden with a set of control pointssatisfying the C^(r) continuity constraints:p=C ^(T) p ^(ƒ),  (19).

In addition, the resulting C^(r) mesh recovers the original NURBSboundary exactly.

Now the C^(r) geometric map G(ξ):{circumflex over (Ω)}→Ω can be obtainedin terms of rational C^(r) basis functions ψ_(i)(ξ), or equivalently,the rational C⁰ Bernstein basis functions ϕ_(j)(ξ) as:

$\begin{matrix}{{{G(\xi)} = {{\sum\limits_{i}^{m}\;{p_{i}^{f}{\psi_{i}(\xi)}}} = {\sum\limits_{j}^{n}{p_{j}{\phi_{j}(\xi)}}}}},} & (20)\end{matrix}$

where m and n are the dimension of the space S_(d) ^(r) and S_(d) ⁰respectively.

Isogeometric Analysis Using rTBS Elements

Following herein, a method of isogeometric analysis using rTBS elementsis described, where the classical Galerkin formulation is applied. Theproblems considered in the following include linear elasticity andPoisson problem.

The governing equation for the linear elasticity is:

$\begin{matrix}\left\{ \begin{matrix}{{{\nabla{\cdot \sigma}} + b} = 0} & {{on}\mspace{14mu}\Omega} \\{\sigma = {D{\nabla_{s}u}}} & \; \\{{\sigma \cdot n} = t} & {{on}\mspace{14mu}\Gamma_{t}} \\{u = \overset{\_}{u}} & {{{on}\mspace{14mu}\Gamma_{u}},}\end{matrix} \right. & (21)\end{matrix}$

where D is the elasticity matrix, b and t refer to body force andtraction respectively, u is the displacement, Γ_(t) and Γ_(u) are theportions of the boundary, where traction and displacement are specified,respectively. The Poisson problem is defined as:

$\begin{matrix}\left\{ \begin{matrix}{{- {\nabla^{2}u}} = f} & {{{in}\mspace{14mu}\Omega},} \\{u = \overset{\_}{u}} & {{{on}\mspace{14mu}\Gamma},}\end{matrix} \right. & (22)\end{matrix}$

where ƒ:Ω→

is a given function and ū denotes prescribed boundary values.

Using the basis constructed in the previous section, we approximate thesolution in the corresponding parametric domain as:

$\begin{matrix}{{{\hat{u}(\xi)} = {{\sum\limits_{i}{u_{i}{\psi_{i}(\xi)}}} = {u^{T}\psi}}};} & (23)\end{matrix}$

where u_(i) corresponds to the approximate solution's Bézier ordinate atthe i-th domain point in the parametric domain {circumflex over(Ω)}_({circumflex over (T)}), as illustrated in FIG. 1. The solutionu(x) over the domain Ω_(T) in the physical space is obtained bycomposing û(ξ) with the inverse of the geometric mapping G⁻¹ such thatu(x):Ω

², where:u(x)={circumflex over (u)}(ξ)∘G ⁻¹(x).  (24).

After inserting the approximate solution and basis functions into thecorresponding weak form of the PDE, we obtain the following mass andstiffness matrices, respectively:M ₀=∫_(Ω) ϕ·ϕdΩ,  (25);K ₀=∫_(Ω) ∇ϕ·∇ϕdΩ,  (26);

for C⁰ elements. For C^(r) elements. The mass and stiffness matrices arecalculated using the fact that the C^(r) basis ψ are linear combinationsof C⁰ basis ϕ:M _(r)=∫_(Ω) ψ·ψdΩ=∫ _(Ω)(Cϕ)·(Cϕ)dΩ=C ^(T) {tilde over (M)} ₀ C,  (27);K _(r)=∫_(Ω) ∇Ω·∇ψdΩ=∫ _(Ω)(C∇ϕ)·(C∇ϕ)dΩ=C ^(T) {tilde over (K)} ₀C,  (28);

where M₀ and {tilde over (K)}₀ are the mass and stiffness matricesrespectively for the same C^(r) elements in terms of the C₀ basis ϕ. Thedifference between {tilde over (M)}₀ and M₀, {tilde over (K)}₀ and K₀ isdue to the potential relocation of the control points to satisfy theC^(r) continuity constraints for the C^(r) elements. The assemblyprocess for such matrices is different from the one described in N.Jaxon, X. Qian, Isogeometric analysis on triangulations, Computer-AidedDesign 46 (2014) 45-57, where the entries in M_(r) and K_(r) arecalculated directly after identifying the basis function ψ_(i)supporting each element from the continuity matrix C. Instead, {tildeover (M)}₀ and {tilde over (K)}₀ are assembled first using the C⁰ basisϕ as in classical C⁰ FEM and then multiplied by the continuity matrix Cto obtain M_(r) and K_(r) as shown in equations (27) and (28). Thisimplementation can be readily applied in any existing FEM routinewithout changing the assembly process. The numerical integration isperformed in each element (micro-element if split is used) by usingstandard and collapsed Gaussian quadrature rules on the boundaries andelement interiors respectively. Specifically, the integrals are firstpulled back onto the parametric domain and then onto a parent element ofright-angled isosceles triangle, as shown in FIG. 1.

Due to the use of boundary MDS mentioned earlier, the Dirichlet boundaryconditions can be imposed similarly as in NURBS based IGA. Typicalstrategies include the least square method (used herein), weakimposition using Lagrange multiplier, and an improved transformationmethod.

Approximation Property of the rTBS Space

It is already known that the set of Bernstein basis functions of degreed over a triangulation form a basis for the polynomial space of degreed. Thus for C⁰ basis, the same error estimate holds as for classicalfinite element methods and optimal convergence rates can be guaranteed.In this section, we focus on how well C^(r) continuous functions can beapproximated by rTBS.

It has been proven in M. Lai, L. Schumaker, Spline functions ontriangulations, Vol. 110, Cambridge University Press, 2007, which isincorporated herein by reference in its entirety, that, if there existsa subspace of S_(d) ^(r)({circumflex over (T)}) with a stable localminimal determining set, then S_(d) ^(r)({circumflex over (T)}) has theapproximation power up to d+1. That is, for every ƒϵH^(d+1), thereexists a spline sϵS_(d) ^(r)({circumflex over (T)}) such that:

$\begin{matrix}{{{f - s}}_{W^{k,{d + 1}}(\overset{\Cap}{\Omega})} \leq {{Ch}_{\overset{\Cap}{T}}^{d + 1 - k}{f}_{W^{{d + 1},{d + 1}}(\overset{\Cap}{\Omega})}\mspace{14mu} 0} \leq k \leq {d.}} & (29)\end{matrix}$

where H^(k) and W^(k,p) are the Hilbert and Sobolev spaces respectivelywith ∥⋅∥ the associated seminorm, h_({circumflex over (T)}) is thelength of the longest edge in T, the constant C depends only on d andthe smallest angle in T, and the Lipschitz constant of the boundary ofΩ.

Although the result relative to equation (29) is derived in the TBSspace, following NURBS based IGA, a similar result in the rTBS space canbe derived as:

$\begin{matrix}{{{{f - {\Pi_{\mathcal{S}}f}}}_{H^{k}{(\tau)}} \leq {C_{w}h_{\overset{\Cap}{T}}^{d + 1 - k}\;{f}_{W^{{d + 1},{d + 1}}(\overset{\Cap}{\Omega})}}},{\forall{f \in H^{d + 1}}},{0 \leq k \leq {d.}}} & (30)\end{matrix}$

where Π_(S) is the projector on the rTBS space S_(d) ^(r), the constantC_(w) differs from C by the extra dependence on the weight function wand τ is an element in the triangulation T. Finally we define theprojector Π_(u):L²(Ω)→u_(d) ^(r) as:Πuƒ:=(Π_(S)(ƒ∘G))∘G ⁻¹ ,∀ƒϵL ²(Ω),  (31);

where u_(d) ^(r) is the space of rTBS on the physical domain Ω (thepush-forward of the rTBS space S_(d) ^(r) on the parametric domain Ω),as shown in FIG. 7.

Now the error estimate on the physical domain Ω can be derived as:

$\begin{matrix}{{{{f - {\Pi_{\mathcal{U}}f}}}_{H^{k}{(T)}} \leq {C_{w}h_{T}^{d + 1 - k}{\sum\limits_{i = 0}^{d + 1}{{{\nabla G}}_{L^{\infty}{({G^{- 1}{(T)}})}}^{i - d - 1}{f}_{H^{i}{(T)}}}}}},{\forall{f \in {{H^{d + 1}(\Omega)}.}}}} & (32)\end{matrix}$

where G is the geometric map, T=G({circumflex over (T)}), and h_(T) isthe longest element edge in T. Thus, equation (32) implies the rTBSspace on the physical domain delivers the optimal rate of convergence,which is d+1−k in terms of the error norm in H^(k) space for polynomialsof degree d, provided that there is a set of local stable basis forS_(d) ^(r) and the geometric map G remains the same for different meshsize h_(T).

Pre-Refinement Smooth Geometric Map

As shown above, in order to evaluate the convergence rate uponh-refinement, the geometric map must remain the same during refinement.As will follow, the need for a smooth pre-refinement geometric map willbe shown. Then a strategy to construct a pre-refinement geometric mapwill be introduced that possesses sufficient smoothness to maintain theconsistency of the geometric map for all subsequent refinements. Thus,the basis for a “smooth-refine-smooth” approach, in accordance with thepresent disclosure, will be provided.

Need for a Pre-Refinement Smooth Geometric Map

Let Δ₀⊂Δ₁⊂Δ₂⊂ ⊂Δ_(n) be a nested sequence of triangulations, whereΔ_(k-1)⊂Δ_(k) means Δ_(k) is a refinement of Δ_(k-1) by subdividing eachtriangle in Δ_(k-1) into several sub-triangles. We denote the C^(r)spline space defined on Δ_(k) as S_(d) ^(r)(Δ_(k)). If the meshsequences in the parametric domain are nested, that is, {circumflex over(T)}₀⊂{circumflex over (T)}₁⊂{circumflex over (T)}₂⊂ . . . ⊂{circumflexover (T)}_(n), and the resulting triangulations in the physical domainunder the geometric map G_(k):{circumflex over (T)}_(k)

T_(k) are also nested as T₀⊂T₁⊂T₂⊂ . . . ⊂T_(n), then the geometric mapG_(k) for the space S_(d) ^(r)(Δ_(k)) is said to be consistent duringthe refinement sequence.

For accurate isogeometric analysis with C^(r) rTBS elements, elementsneed to be sufficiently small. One way is to first obtain C⁰ coarse meshfrom the procedure outlined above with respect to the construction of aC⁰ geometric map G₀, perform uniform refinement on the C⁰ mesh until theelements are small enough for accurate analysis, and then impose C^(r)smoothness constraints through a macro-element based DC technique or GE.That is, for space S_(d) ^(r), first create a nested sequence oftriangulations T₀⊂T₁⊂T₂⊂ . . . ⊂T_(n) in the physical domain, and thenimpose C^(r) continuity constraints on each mesh T_(k) and relocate thedependent control points to ensure C^(r) smoothness. Although such a“refine-then-smooth” approach is able to create arbitrarily small rTBSelements with desired continuity for analysis, the resulting geometricmap may not be consistent during the refinement. That is, the resultingelements are not nested after the refinement. As such, optimalconvergence cannot be achieved with such an approach. To demonstratesuch a lack of consistency in the resulting C^(r) geometric map with the“refine-then-smooth” approach, two examples will be discussed.

Lack of Consistency Example 1

The notation in various refinement and splits are as follows. Thesubscripts u,ct,ps are used to indicate uniform refinement, CT split andPS split, respectively, which are used in the same sequential order asthese refinements are performed. The superscript indicates the order ofsmoothness. For example, T_(ct,u,u,ps) ¹ represents a C¹ smooth meshobtained by performing a CT split followed by two uniform refinementsand a PS split on the mesh T. Note that the CT and PS split are usuallyfollowed by imposing continuity constraints. Here ct* and ps* are usedto indicate the respective split without imposing continuityconstraints.

An example is given below where a domain is initially parameterized intofive C⁰ elements, as shown in FIGS. 8A and 8B. Also, FIGS. 9A-9F shows anested sequence of triangulations from such initial C⁰ elements throughuniform refinement. CT splits are then performed to obtain splines inthe S₃ ¹ FIGS. 9A-9C are C⁰ meshes before continuity constraints areactually imposed, where the free control points 90 are points that canbe chosen freely, but continuity-satisfied, dependent control points 92and non-continuity-satisfied, control points 94 are points that aredetermined through the continuity constraints. That is,continuity-satisfied, dependent control points 92 have already satisfiedthe continuity condition, while the non-continuity-satisfied, dependentcontrol points 94 are not and need to be relocated to obtain C¹smoothness. As can be seen, the non-continuity-satisfied, dependentcontrol points 94 are distributed only near the common edges shared bythe initial five elements in FIG. 8B. After relocating the controlpoints, as shown in FIGS. 9D-9F, the meshes have been locally changedand are no longer nested. The relocation of some control points p_(i)(the non-continuity-satisfied, dependent control points 94) to satisfythe continuity constraints, thus, leads to a change of the geometric mapaccording to equation (20).

Lack of Consistency Example 2

The second example concerns a kind of macro-elements in the supersplinespace S_(d) ^(r,ρ), ρ>r, as defined in equation (11), wheresuper-smoothness C^(ρ) happens at the vertices or edges of themacro-triangles. Uniform refinement of such elements followed by thesame macro-element technique to achieve super-smoothness atmacro-element vertices or edges would lead to inconsistent geometricmaps. An example is given in FIGS. 10A-10C. Specifically, the initialmesh is illustrated in FIG. 10A in superspline space S₅ ^(1,2) obtainedby DC based on the quintic C¹ macro-element technique. This initial meshis globally C¹ but with C² smoothness at the vertices ofmacro-triangles. If this initial T^(1,2) mesh is refined as shown inFIG. 10B where the smoothness at the points 100 are C² and thesmoothness at points 102 in the refined mesh is still C¹. In order toobtain stable basis in the S₅ ^(1,2) with the same macro-elementtechnique, all macro vertices need to be C² smooth. Thus, the controlpoints 104 need to be relocated to achieve C² smoothness at all macrovertices, since many of the control points only possess C¹ smoothnessbefore refinement. As the mesh is refined, as illustrated in FIG. 10C,there are more and more vertices that need to be relocated to achieve C²smoothness. Such relocation of control points to achieve higher-ordersmoothness at more and more vertices and edges of refinedmacro-triangles, due to super-smoothness requirement at vertices ofmacro-triangles in the S₅ ^(1,2) macro-element space, leads to the lossof consistency of the geometric map.

Solution

To overcome such inconsistency of the geometric map during therefinement, a geometric map with sufficient smoothness may beconstructed before the refinement. For the usual C^(r) splines, thepre-refinement geometric map may be at least C^(r) smooth. Forsuperspline space S_(d) ^(r,ρ) where super-smoothness occurs at thevertices or edges of macro-triangles, a C^(ρ) pre-refinement map may beconstructed before the refinement. In this way, all subsequent refinedelements are globally C^(ρ) smooth and the super-smoothness required atthose vertices and edges are therefore satisfied. The refinementsequence is nested and the geometric map remains unchanged.

Note if super-smoothness happens at the interior vertices or edges ofmacro-triangles after the splits, refinement of such elements does notneed to relocate dependent control points to satisfy the continuityconstraints. This is because, at these internal vertices and edges of amacro-triangle, the continuity is already C^(∞). For example, althoughthe quadratic C¹ PS macro-element and cubic C¹ CT macro-element spacesare also superspline spaces, their super-smoothness occurs inside themacro-elements where the geometric smoothness is infinity. Thus, theycan be treated as regular C^(r) spaces in terms of smoothnessrequirement of the pre-refinement geometric map.

Construction of a Pre-Refinement Smooth Geometric Map

To control or avoid the inconsistency of the C^(r) geometric map duringthe refinement, a pre-refinement map may be constructed. Thispre-refinement map may be designed to meet a sufficiency criteria withrespect to smoothness. The pre-refinement map can have refinementsperformed on it, and the C^(r) continuity constraints can be followed toobtain stable basis for C^(r) rTBS elements, such as discussed withrespect to construction of C^(r) spline basis ψ(ξ) on {circumflex over(T)}. The refined mesh inherits the continuity. The refined controlpoints, therefore, do not need to be relocated since they alreadysatisfy the continuity conditions. Thus, the resulting meshes are C^(r)smooth and nested, and the geometric map remains the same for allsubsequent refinements.

A sufficiently smooth pre-refinement mesh for a given domain can beobtained, such as described above. For example, if the domain is boundedby straight line segments and recalling that the domain points in theparametric mesh satisfies the smoothness condition, the control pointscan be generated in the location obtained by an affine transformation ofthe domain points in the parametric mesh. In this case, the physicalmesh obtained will satisfy the needed smoothness condition as well. Ifthe domain has curved boundaries, a C^(r) mesh can be constructed with aset of stable local basis (such as described above with respect toconstruction of C^(r) spline basis ψ(ξ) on or the control points can berelocated to satisfy the smooth conditions by Gaussian elimination.

It is worthy to mention the difference between the pre-refinement C^(r)smooth geometric map described here and the C^(r) map described abovewith respect to construction of C^(r) spline basis ψ(ξ) on T. The reasonfor a pre-refinement smooth geometric map is geometric. That is, a C^(r)map is created that maps the parametric domain to the physical domainand recovers the original boundary, without the need for a set of stablelocal basis for analysis. Such a smooth pre-refinement map is designedso that, during the refinement, the control points do not need torelocate since they would already satisfy the continuity conditions.Thus, the pre-refinement map is kept consistent during the refinement.On the other hand, in order to approximate a field in analysis, a set ofstable basis should be constructed by the methods discussed with respectto construction of C^(r) spline basis ψ(ξ) on {circumflex over (T)}. Forexample, in cubic space S₃ ¹(T_(ct)) with CT split, a C¹ pre-refinementmap can be created using CT split, as shown in FIG. 11A. Even if it isC¹ after uniform refinement, as illustrated in FIG. 11B, the CT split isperformed again in order to construct a set of stable local basis in S₃¹(T_(ct)) to be used in analysis, as shown in FIG. 11C. Thus, FIGS.11A-11C illustrate smooth-refine-smooth procedure, where the purpose ofthe first CT split is to construct a pre-refinement smooth geometricmap, while the CT split at the last step is to construct a set of C¹stable local basis for analysis.

Numerical Results

In the following discussion, “optimal convergence” rates are discussed,as is achieving “optimal convergence” using rTBS-based isogeometricanalysis for different problems with different elements. First, thiswill be illustrated with respect to a domain bounded by straight linesegments, where optimal convergence rates are achieved in all presentedspaces. Then, with respect to domains with curved boundaries, if nopre-refinement smooth map is constructed to keep the C^(r) geometric mapconsistent during refinement, the convergence rates will be shown to belower than optimal. Then, with the pre-refinement smooth map asdiscussed above with respect to pre-refinement smooth geometric map,optimal convergence rates will be discussed. The results providedhereafter are demonstrated on two examples. A first is a Poisson problemon a complex domain with three holes and one elastic problem on a platedomain. We also show the advantage of local refinement in rTBS for theelasticity problem. Note, in all examples below, the element size in theconvergence study refers to the maximal length of the edges in themicro-elements.

Domain with all Straight Boundaries: Poisson Problem

The first example is a triangular domain with a triangular hole as shownin FIGS. 12A-12F. All sides consist of straight line segments. Thegoverning equation is equation (22) with the open set domain Ω beingdefined as:

$\begin{matrix}{\Omega:=\left\{ \left( {x,y} \right) \middle| {\left\lbrack {{{{{\left( {0 < x \leq \frac{3}{2}} \right)\&}\left( {0 < y < {3\sqrt{3x}}} \right)}\bigcup\left( {\frac{3}{2} < x < 1} \right)}\&}\left( {0 < y < {3\sqrt{3}\left( {1 - x} \right)}} \right)} \right\rbrack\backslash\;{\left. \quad\left\lbrack {{{{{\left( {\frac{9}{8} \leq x \leq \frac{3}{2}} \right)\&}\left( {\frac{3\sqrt{3}}{8} \leq y \leq {3\sqrt{3}\left( {x - \frac{1}{4}} \right)}} \right)}\bigcup\left( {\frac{3}{2} < x \leq \frac{15}{8}} \right)}\&}\left( {\frac{3\sqrt{3}}{8} \leq y \leq {3\sqrt{3}\left( {\frac{3}{4} - x} \right)}} \right)} \right\rbrack \right\}.}} \right.} & (33)\end{matrix}$

The body force is:ƒ(x,y)=−(x ² +y ²)e ^(xy),  (34).

and the exact solution is given by:u(x,y)=e ^(xy),  (35).

Based on our parameterization strategy in described above, theparametric domain illustrated in FIG. 12A is the same as the physicaldomain, since the boundary edges are all straight. The initial linearmesh is obtained by the Delaunay triangulation of the problem domain, asshown in FIGS. 12B and 12C. By degree elevation, meshes corresponding toquadratic, cubic and quintic C⁰ spline spaces are obtained, as shown inFIGS. 11C, 11D, and 11E, respectively. Only the physical meshes areshown since the parametric meshes are the same as the correspondingphysical meshes.

Meshes corresponding to C^(r) spaces with stable basis for analysis canbe obtained through the macro-element techniques via either the DC or GEmethod as described above. For example, FIGS. 13A-13H show the C^(r)mesh corresponding to C⁰ meshes discussed above with respect to FIGS.12A-12F, including quadratic C¹ PS macro-element space S₂ ¹(T_(ps)),cubic C¹ CT macro-element space S₃ ¹(T_(ct)), quintic C¹ polynomialmacro-element space S₅ ¹(T), S₅ ^(1,2)(T) and quintic C² PSmacro-element space S₅ ²(T_(ps)) and S₅ ^(2,3)(T_(ps)). In FIGS.13A-13H, solid dots represent free nodes and hollow dots representdependent nodes. The dimension of each space (number of independentbasis functions) is also reported. As can be seen, for S₂ ¹(T_(ps)) andS₃ ¹(T_(ct)), although the MDS obtained by DC and GE may be different,the dimension of the spaces are exactly the same, and so are thenumerical solutions resulted from the two methods. For the quintic C¹and C² spaces, DC yields superspline spaces which have smaller dimensionthan the ones obtained by GE. For example, the space S₅ ^(1,2)(T)obtained by DC has dimension 96 (FIG. 13E) while the space S₅ ¹(T)obtained by GE has dimension 120 (FIG. 13F), and the space S₅^(2,3)(T_(ps)) obtained by DC has dimension 120 (FIG. 13G) while thespace S₅ ²(T_(ps)) obtained by GE has dimension 252 (FIG. 13H).

To study the convergence, uniform refinements are performed on theinitial C⁰ meshes (FIGS. 12D-12F) before the same macro-elementtechniques are used to obtain stable C^(r) basis. The refinementsequences for these spaces are {T_(ps),T_(u,ps),T_(u, . . . ,u,ps), . .. } in S₂ ¹(T_(ps)), in {T_(ct),T_(u,ct),T_(u, . . . ,u,ct), . . . } inS₃ ¹(T_(ct)), {T,T_(u),T_(u, . . . ,u,) . . . } in S₅ ^(1,2)(T),{T,T_(u),T_(u, . . . ,u), . . . } in S₅ ¹(T_(ps)),{T_(ps),T_(u,ps),T_(u, . . . ,u,ps), . . . } in S₅ ^(2,3)(T_(ps)) and{T_(ps),T_(u,ps),T_(u, . . . ,u,ps), . . . } in S₅ ²(T_(ps)). Note, inthis example, no pre-refinement smooth map was explicitly constructed.The reason for this is that, for domains bounded by straight linesegments, the parametric mesh are identical to the physical mesh and thegeometric map is in fact C^(∞) smooth.

In the convergence study, we compute the analysis error by:

$\begin{matrix}{{e_{u} = \left\lbrack {\int_{\Omega}{{\left( {u_{num} - u_{exact}} \right) \cdot \left( {u_{num} - u_{exact}} \right)}{\mathbb{d}\Omega}}} \right\rbrack^{1/2}},} & (36)\end{matrix}$

where u_(num) and u_(exact) are the numerical and exact solutions,respectively. The longest edge h_(max) of the triangles in the physicalmesh is considered as the mesh parameter. FIGS. 14A-14C illustrate theerror measured in the L²-norm versus mesh parameter. As shown in FIGS.14A-14C, optimal convergence rates are achieved in the tested C⁰, C¹ andC² spaces, including in the superspline spaces S₅ ^(1,2)(T) in FIG. 14Band S₅ ^(2,3)(T_(ps)) in FIG. 14C. Particularly, comparing with theregular C^(r) spline space, the same optimal convergence rate isobtained in the superspline space, but with far fewer degree offreedoms. For example, for space S₅ ²(T_(ps)), the number of free DOFsobtained via the GE method are respectively 252, 843, 3063, and 11679,corresponding to the refinement sequences in FIG. 13C. On the otherhand, for the superspline S₅ ^(2,3)(T_(ps)), the number of DOFs obtainedvia the DC method are respectively 120, 360, 1200, 4320 and 16320,corresponding to the refinement sequences in FIG. 13C. Note that forspaces S₂ ¹(T_(ps)) and S₃ ¹(T_(ct)), DC and GE yields exactly the sameanalysis results and their convergence curves overlap with each other inthis figure.

The refinement steps, the methods for obtaining C^(r) stable basis foranalysis, and the convergence rates for each type of elements aresummarized in Table 1. It can be seen that optimal convergence rateshave been achieved with all types of C^(r) rTBS elements.

TABLE 1 smooth-refinement-smooth steps and convergence rates ofdifferent C^(r) spaces for the problem in FIG. 12A-12F. StepPre-refinement map Refine- Stable basis Conv. Space Smoothness SplitMethod ment Split Method rate S₂ ¹(T_(ps)) C^(∞) uniform PS DC or 3.0 GES₃ ¹(T_(ct)) C^(∞) uniform CT DC or 3.8 GE S₅ ¹(T) C^(∞) uniform GE 6.0S₅ ^(1,2)(T) C^(∞) uniform DC 6.0 S₅ ²(T_(ps)) C^(∞) uniform PS GE 5.8S₅ ^(2,3)(T_(ps)) C^(∞) uniform PS DC 5.9

Domain with Curved Boundaries: Poisson Problem

In this example we solve a Poisson problem on a L-shaped domain withthree holes, as shown in FIG. 15A. White nodes represent the controlpoints. The boundaries of the domain are represented in NURBS withweights so chosen that exact circular holes are represented. Thegoverning equation is equation (22) with the open set domain beingdefined as:Ω:={((x,y)|[0≤x≤16)&(0≤y≤1.6)]\[((8<x<16)&(8<y<16))∪((x−4)²+(y−4)²<4)∪((x−12)²+(y−4)²<4)∪((x−4)²+(y−12)²<4)]}.  (37).

The body force is:ƒ(x,y)=2 sin(x)sin(y),  (38);

and the exact solution is given by:u(x,y)=sin(x)sin(y).  (39).

To create the initial parametric mesh, we first extract quadratic Béziercurves from the NURBS boundary curves of the physical domain.Particularly, four rational Bézier segments are exacted from eachcircular boundary. The end points of these Bézier curves are connectedto form the initial parametric domain, which is then triangulated toobtain the initial parametric mesh, as shown in FIG. 15B. We thenreplace the boundary control points of the parametric mesh withcorresponding points from the physical boundary to obtain the initialphysical mesh, as shown in FIG. 15C. The C^(r) stable basis for analysiscan be obtained as usual by degree elevation followed by either the DCor GE method based on the macro-element techniques. For example, FIGS.16A-16H show parametric and physical meshes for different C^(r) spaces.In particular, the meshes correspond to quadratic C¹ PS macro-elementspace S₂ ¹(T_(ps)), cubic C¹ CT macro-element space S₃ ¹(T_(ct)),quintic C¹ polynomial macro-element space S₅ ¹(T) and S₅ ^(1,2)(T).

The refinement sequences used to evaluate the convergence are{T_(ps),T_(ps,u,ps),T_(ps,u, . . . ,u,ps), . . . } in S₂ ¹(T_(ps)),{T_(ct),T_(ct,u,ct),T_(ct,u, . . . ,u,ct), . . . } in S₃ ¹(T_(ct)),{T,T_(u),T_(u, . . . ,u), . . . } in S₅ ¹(T), and{T_(ps),T_(ps,u),T_(ps,u, . . . ,u), . . . } in S₅ ^(1,2)(T). Themethods for constructing pre-refinement smooth maps, refinement, andmethods of basis construction in each C^(r) space are summarized inTable 2. For example, for the superspline space S₅ ^(1,2)(T), theinitial physical mesh is C² smooth obtained by GE with PSmacro-elements. The four pairs of figures (i.e., FIGS. 16/16B, 16C/16D,16E/16F, 16G/16H) correspond to the four initial pairs of parametric andphysical meshes in the four convergence studies, respectively. Theconvergence rates obtained in C⁰ and C¹ spaces are all optimal as shownin FIGS. 17A and 17B, where the convergence rates for quadratic, cubicand quintic elements are 3, 4, and 6 respectively in all C⁰, C¹ andC^(1,2) spaces.

TABLE 2 Smooth-refinement-smooth steps and convergence rates ofdifferent C^(r) spaces for the problem in FIG. 15A-C. StepPre-refinement map Refine- Stable basis Conv. Space Smoothness SplitMethod ment Split Method rate S₂ ¹(T_(ps)) C¹ PS DC uniform PS DC 3.1 S₃¹(T_(ct)) C¹ CT DC uniform CT DC 3.9 S₅ ¹(T) C¹ GE uniform GE 6.0 S₅^(1,2)(T) C² PS GE uniform DC 5.7

If the refine-then-smooth strategy is used and, thus, no smoothpre-refinement map is constructed, the geometric map changes andconvergence rate decreases as the mesh is refined. As shown in FIG. 17C,although the rate in S₂ ¹(T_(ps)) is optimal, the rates in S₃ ¹(T_(ct))and S₅ ¹(T) are 3, and the rate in S₅ ^(1,2)(T) decreases quickly toabout 4.1. Note that the extreme large errors in S₅ ^(1,2)(T) on for thecoarse meshes are due to the poor mesh quality.

Domain with Curved Boundaries: Linear Elasticity

In the third example, we apply our approach to a well-known linearelasticity problem: an infinite plate with a circular hole underconstant in-plane tension in the x-direction. The infinite plate ismodeled by a finite quarter plate as shown in FIG. 18 with the governingequation (21). L is the length of the edge, R is the radius of thecircle and τ is the thickness of the plate. E and ν represent theYoung's modulus and Poisson ratio, respectively. The exact solution,evaluated at the boundary of the finite quarter plate, is applied as aNeumann boundary condition.

The initial parameterization of the physical domain is shown in FIGS.19A-19C. The NURBS boundary curves of the given domain are firstextracted as rational quadratic Bézier curves. After connecting the endpoints of the Bézier curves, we obtain an initial parametric domain andtriangulate it, as shown in FIG. 19B. Then, we replace the boundarycontrol points with corresponding points on the physical boundary toobtain the initial physical mesh as shown in FIG. 19A. To improve themesh quality and analysis results, we use the same smoothed parametricmesh, as in FIG. 19C, by minimizing the difference of the internalcorner angles between the parametric and physical domains.

Through degree elevation on the parametric and physical meshes in FIGS.19A-19C, we obtain cubic and quintic C⁰ meshes in S₃ ⁰(T) and S₅ ⁰(T),as shown in FIG. 20A-20P. Meshes in C^(r) spaces for analysis are alsoobtained by either DC or GE method, as shown in FIGS. 20A-20P, wheresolid dots represent free nodes and hollow dots represent dependentnodes.

The refinement sequences used to evaluate the convergence are{T,T_(u),T_(u, . . . ,u), . . . } in all C⁰ spaces,{T_(ps),T_(ps,u,ps),T_(ps,u, . . . ,u,ps), . . . } in S₂ ¹(T_(ps)),{T_(ct),T_(ct,u,ct),T_(ct,u, . . . ,u,ct), . . . } in S₃ ¹(T_(ct)),{T,T_(u),T_(u, . . . ,u), . . . } in S₅ ^(1,2)(T),{T,T_(u),T_(u, . . . ,u), . . . } in S₅ ¹(T),{T_(ps),T_(u,ps),T_(u, . . . ,u,ps), . . . } in S₅ ^(2,3)(T_(ps)) and S₅²(T_(ps)). The steps of convergence analysis and methods of basisconstruction in each C^(r) space is summarized in Table 3. Pairsrepresented by FIGS. 20A/20B, 20C/20D, 20E/20F, and 20G/20H correspondto the four initial pairs of parametric and physical meshes in the fourconvergence studies respectively. Note that the pre-refinement smoothgeometric map, C² and C³ used in space S₅ ²(T_(ps)) and S₅^(2,3)(T_(ps)) are both obtained by GE after imposing continuityconstraints on the quintic C⁰ mesh.

TABLE 3 Smooth-refinement-smooth steps and convergence rates ofdifferent C^(r) spaces for the problem in FIG. 18. Step Pre-refinementmap Refine- Stable basis Conv. Space Smoothness Split Method ment SplitMethod rate S₂ ¹(T_(ps)) C¹ PS DC uniform PS DC 2.0 S₃ ¹(T_(ct)) C¹ CTDC uniform CT DC 2.8 S₅ ¹(T) C¹ GE uniform GE 5.0 S₅ ^(1,2)(T) C² GEuniform DC 5.0 S₅ ²(T_(ps)) C² GE uniform PS GE 5.0 S₅ ^(2,3)(T_(ps)) C³GE uniform PS DC 4.7

The energy error is evaluated by:

$\begin{matrix}{{e_{stress} = \left\lbrack {\frac{1}{2}{\int_{\Omega}{{\left( {ɛ_{num} - ɛ_{exact}} \right) \cdot D \cdot \left( {ɛ_{num} - ɛ_{exact}} \right)}{\mathbb{d}\Omega}}}} \right\rbrack^{1/2}},} & (40)\end{matrix}$

where ε_(num) and ε_(exact) are the numerical and exact strain vectorsrespectively. The mesh parameter is evaluated as the longest edgeh_(max) of the triangles in the physical mesh. Again, optimalconvergence rates are achieved in C⁰, C¹ and C² spaces, as shown inFIGS. 21A-21C where quadratic, cubic, and quintic rates are obtained forthe energy norm error using quadratic, cubic, and quintic elements,respectively. In the superspline spaces S₅ ^(1,2)(T) and S₅^(2,3)(T_(ps)), optimal rates are also observed as shown in FIGS. 21Band 21C. This demonstrates the efficiency of supersplines for analysison a per-node basis since far fewer degree of freedoms are used insuperspline space than the regular spline space of the same degree.Also, FIG. 22 provides a plot of the convergence curve on a node basis.It can be seen that at the same number of nodes in three quinticelements, S₅ ⁰, S₅ ¹, S₅ ^(2,3), higher continuity leads to smallererrors.

However, if we use the refine-then-smooth strategy, that is, withoutconstructing the smooth pre-refinement map, the local modification ofthe control points to obtain the C^(r) basis would change the geometricmap. Consequently, the convergence rates are reduced, as shown in FIG.23. The convergence rates in all C¹ and C² spaces with polynomialdegrees ranging from quadratic to quintic are about the same, whichrange from 1.4 to 1.6 and are far below the optimal values.

One advantage of using triangular meshes in analysis is the ease oflocal refinement of meshes. We use the Rivara method, which uses theelement-wise strain energy error to guide the local refinement. Theelements with large error are bisected across one of their edges. Thelocal refinement of S₃ ¹(T_(ct)) mesh is shown in FIGS. 24A-24D. Thecomparison of error measured in the L²-norm of stress vs. number ofdegree of freedoms between uniform and adaptive refinement is shown inFIG. 25. The local refinement sequence shows superior advantage withsame error obtained using only one third of degree of freedoms ofuniform refinement. Clearly, the local refinement exhibits superioradvantage over uniform refinement by leading to the same accurateresults with fewer number of degree of freedoms. Specifically, to obtainthe same order of error, the uniform refinement requires about threetimes the number of degree of freedoms as much as the adaptiverefinement.

As described above, the smooth-refine-smooth approach to rTBS basedisogeometric analysis can achieve convergence rates for all C^(r) rTBSelements. This convergence can be considered “optimal” relative to agiven criteria. For example, “optimal convergence” may be defined, asdescribed above, where equation (32) implies the rTBS space on thephysical domain delivers the optimal rate of convergence, which is d+1−kin terms of the error norm in H^(k) space for polynomials of degree d,provided that there is a set of local stable basis for S_(d) ^(r) andthe geometric map G remains the same for different mesh size hr.Likewise, “sufficient smoothness” may be defined, for example, as thesmoothness needed to maintain the consistency of the geometric map forsubsequent refinements.

As described above, for a given NURBS bounded domain with arbitrarytopology, the rTBS based parameterization can be fully automated, suchthat user intervention is not required to select or complete refinement.Various sets of globally C^(r) continuous basis can be constructed byimposing continuity constraints on adjoining triangle elements, througheither DC or GE method. Error estimates indicate that the constructedC^(r) space delivers optimal convergence rates, provided the geometricmap remains the same during refinement.

Referring to FIG. 26, a computer network 1000 is illustrated, which maybe configured to carry out a process in accordance with the presentdisclosure. The computer network 1000 may be formed by a plurality ofworkstations 1002 that are connected to one or more servers 1004.Regardless of whether the computer system is a server 1004 orworkstation 1002, an underlying hardware architecture 1006 of thecomputer systems is illustrated. The hardware architecture 1006 mayinclude one or more CPUs 1008, which may include one or more caches1010. The CPU 1008 is generally connected through a bridge 1012 tomemory 1014 and, in some cases, an additional non-local cache 1016. Somesystem may include dedicated graphics processing units (GPUs) 1018 thathave been adapted from processors utilized to simply drive a display1020 to a secondary, specialized processor that the CPU 1008 can utilizeto offload tasks fitting the specialized capabilities of the CPU 1008,such as transcoding operations and many others. In any case, the generalcomputer architecture 1006, regardless of workstation 1002 or server1004 provides a CPU 1008 and memory 1014 and may be supplemented bysecondary processing and memory components, such as a GPU 1018 andvarious caches 1010, 1016 dedicated to particular situations. In thisregard, the above-described components may be conceptualized as aCPU/memory sub-system 1022.

The computer architecture 1006 also includes a bus or multiple buses1024 that connect the above-described CPU/memory sub-system 1022 toother, slower components of the computer architecture 1006. For example,the buses 1024 may provide connections to a universal serial bus (USB)hub or controller 1026 and/or dedicated, bus-connected I/O devices 1028.Of course, I/O connections may vary substantially; however, in allcases, the bus 1024 provides connections to one or more hard drives1030. These hard drives 1030 may take many forms and, more recently,include hardware advances such as solid-state drives, but are uniformlypresent in workstations or personal computers 1002 and servers 1004.This is because traditional notions of computer architecture can beconceptualized as, at a minimum, including a CPU/memory sub-system 1022and a mass-storage sub-system 1032.

Thus, the above-described computer systems may be used to implement theabove-described techniques. For example, referring to FIG. 27, aflowchart is provided that sets forth non-limiting, examples of stepsthat may be carried out, for example, using the system of FIG. 26. Theprocess 1050 begins at process block 1052 by analyzing a CAD model togenerate a pre-refinement geometric map. As described, thepre-refinement geometric map serve as a map between a parametric meshand a physical mesh of the computer aided design model. As such, thepre-refinement geometric map is designed to have a smoothness projectedto maintain a consistency of a mesh based on the pre-refinementgeometric map during a refinement of the mesh. That is, thepre-refinement geometric map has smoothness projected to maintainconsistency of the geometric map between the parametric mesh and thephysical mesh during refinement by making the control points satisfyingthe continuity constraints before refinement so that they do not need tobe relocated after refinement to satisfy the continuity constraints.Thus, the process begins with a “smooth” step.

At process bock 1054, the mesh is refined based on the pre-refinementgeometric map to converge toward a refinement criteria associated withthe CAD model. More particularly, the parametric mesh and the physicalmesh are refined based on the pre-refinement geometric map to convergetoward a refinement criteria associated with the CAD model. Thus, a“refine” step is performed.

At process block 1056, macro-element techniques are used to obtainstable C^(r) smooth basis for analysis. To this end, asmooth-refine-smooth approach is provided. In the above-describedsmooth-refine-smooth approach, the smoothness in the pre-refinementgeometric map keeps the geometric map consistent during the refinement.The smoothness in the last step is used to obtain stable C^(r) basis foranalysis.

Therefore, in order to overcome the inconsistency of geometric mapduring the refine-then-smooth approach, the present disclosure providessystems and methods to construct a pre-refinement map that possessessufficient continuity for subsequent refinements. Thus, the relocationof control points is avoided and the map stays unchanged duringrefinement. By constructing such a pre-refinement geometric map withsufficient smoothness—convergences that meet a given threshold—have beenachieved. In one example a smoothness conditions or threshold forevaluation for the pre-refinement geometric map may be one that is C^(r)smooth for regular C^(r) elements and C^(ρ) smooth in cases ofsuperspline spaces S_(d) ^(r,ρ), ρ>r, where super-smoothness occurs atthe vertices or edges of macro-triangles.

Numerical results verified that convergence rates are improved to beoptimal in different spaces with the introduction of such smoothpre-refinement maps. This demonstrates that C^(r) rTBS elements possesssuperior efficiency on a per-node basis over C⁰ elements. Such nodalefficiency is especially pronounced in the case of supersplines.

It is understood that the present invention is not limited to thespecific applications and embodiments illustrated and described herein,but embraces such modified forms thereof as come within the scope of thefollowing claims.

What is claimed is:
 1. A system for creating a mesh from a computeraided design model during an isogeometric analysis process, the systemcomprising: a memory having access to a computer aided design (CAD)model of an object, wherein the geometry of the CAD model of the objectis described using at least one C^(r) smooth non-uniform RationalB-splines (NURBS) representation, where r≥1; a processor configured tocarry out an isogeometric analysis process that includes: accessing theCAD model of the object from the memory; generating, using the NURBSrepresentation, a pre-refinement geometric map of the CAD object thathas a smoothness projected to maintain a consistency of a mesh based onthe pre-refinement geometric map during a refinement of the mesh suchthat the pre-refinement geometric map is at least C^(r) smooth; defininga plurality of elements by refining the mesh based on the pre-refinementgeometric map to converge toward a refinement criteria associated withthe CAD model; and carrying out the isogeometric analysis of the objectbased on at least the plurality of elements.
 2. The system of claim 1wherein the processor is further configured to further define theplurality of elements by smoothing the mesh after refining to obtain astable C^(r) basis for the isogeometric analysis.
 3. The system of claim2 wherein the processor is further configured to further define theplurality of elements by utilizing a macro-element technique to obtainthe stable C^(r) basis for the isogeometric analysis.
 4. The system ofclaim 1 wherein the processor is further configured to generate thepre-refinement geometric map to serve as a map between a parametric meshand a physical mesh of the computer aided design model.
 5. The system ofclaim 4 wherein the processor is further configured to generate thepre-refinement geometric map to include control points corresponding tofree domain points chosen as free control points.
 6. The system of claim4 wherein the processor is further configured to select the smoothnessprojected to maintain the consistency of the mesh by determining controlpoints that do not need to be relocated during refinement of the meshbecause the control points satisfy a continuity constraint.
 7. Thesystem of claim 6 wherein the processor is further configured to use thecontinuity constraint to obtain a C^(r) basis on domain points in theparametric mesh.
 8. The system of claim 6 wherein the processor isfurther configured to generate the pre-refinement geometric map toinclude dependent control points that satisfy the continuity constraint.9. The system of claim 1 wherein the processor is further configured toselect the smoothness projected to maintain the consistency of the meshby determining control points that do not need to be relocated duringrefinement of the mesh because the control points satisfy a continuityconstraint.
 10. The system of claim 1 wherein the processor is furtherconfigured to evaluate a smoothness condition for the pre-refinementgeometric map including one of determining that the pre-refinementgeometric map is C^(r) smooth for regular C^(r) elements and C^(ρ)smooth in cases of superspline spaces S_(d) ^(r,ρ), ρ>r, wheresupersmoothness occurs at vertices or edges of macro-triangles.
 11. Amethod for creating a mesh from a computer aided design model during aisogeometric analysis process, the method comprising: accessing to acomputer aided design (CAD) model of an object, wherein the geometry ofthe CAD model of the object is described using a t least one C^(r)non-uniform Rational B-splines (NURBS) representation, where r≥1;generating, using the NURBS representation, a pre-refinement geometricmap to serve as a map between a parametric mesh and a physical mesh ofthe computer aided design model, wherein the pre-refinement geometricmap has smoothness projected to maintain consistency of the parametricmesh or the physical mesh during refinement by determining controlpoints that do not need to be relocated during refinement of theparametric mesh or the physical mesh because the control points satisfya continuity constraint, wherein the pre-refinement geometric map is atleast C^(r) smooth; defining a plurality of elements by refining theparametric mesh or the physical mesh based on the pre-refinementgeometric map to converge toward a refinement criteria associated withthe CAD model; and carrying out the isogeometric analysis of the objectbased on at least the plurality of elements.
 12. The method of claim 11further comprising further defining the plurality of elements bysmoothing the mesh after refining to obtain a stable C^(r) basis for theisogeometric analysis.
 13. The method of claim 12 further comprisingfurther defining the plurality of elements by utilizing a macro-elementtechnique to obtain the stable C^(r) basis for the isogeometricanalysis.
 14. The method of claim 12 further comprising generating thepre-refinement geometric map to include control points corresponding tofree domain points chosen as free control points.
 15. The method ofclaim 11 further comprising using the continuity constraint to obtain aC^(r) basis on domain points in the parametric mesh.
 16. The method ofclaim 11 further comprising evaluating an optimal criteria forconvergence.
 17. The method of claim 16 wherein the optimal criteriaincludes evaluating rational triangular Bézier splines (rTBS) space on aphysical domain to deliver an objective rate of convergence.
 18. Themethod of claim 11 wherein the isogeometric analysis process isperformed automatically by a computer system.
 19. The method of claim10, wherein generating the pre-refinement geometric map comprises:subdividing each of a plurality of NURBS curves along a boundary of theobject from the NURBS representation into a set of Bézier curves viaknot insertion; connecting end points of the Bézier curves to form apolygonal parametric domain; and triangulating the polygonal parametricdomain.